It’s been a long time, but we are back again. We previously had discussed how to generate polygons by tracing a circle around a given center and placing vertices at randomly spaced angles and radii. Furthermore, we explained how to identify convexity of a given polygon. Continuing on this thread, we will explore how to determine whether a given point in space is inside or outside a given polygon. For this post, we will look at only convex polygons, but in a future post we will provide code to determine inside/outside for all cases.

First, let’s initialize the function ‘polygonContainsPoint’. The input arguments will be px and py, which are vectors of the x and y coordinates for the polygon. Like our previous examples, the polygon will not be closed, i.e. each vertex is unique (the final vertex is not identical to the first) and the vertices traverse the polygon in order, either counterclockwise or clockwise. The point argument will be a 2-item vector of the (x,y) coordinates for the point in question. The output argument “insidePoly” will be true if the point is within the polygon and false if outside the polygon.

function insidePoly = checkPointToPolygon(px, py, point)

% initialize some variables
insidePoly = true;
numNodes = numel(px);

% Error handling
if numNodes ~= numel(py) || numNodes < 3
    error('Polygon vertices not properly defined.')
end
if numel(point) ~= 2
    error('The point is not defined in 2D space.');
end

In the current code we are discussing, we will only investigate points which are within a convex polygon, therefore we will utilize the previous checkConvex code and run only on appropriate polygons:

% Perform a simple test of inside/outside which applies for convex polygons only
if ~ checkConvex(px, py)
    error('Polygon is concave. Please come back for generic case of algorithm.');
    % TODO: Generic case for both convex and concave 
else
    % Next code block for convex case 
end

For convex polygons, for a point to be inside, it must lie on the same side of each segment of the polygon. Therefore one algorithm is to check each segment in the polygon to see what angle is formed by the point and the segment. One way of doing this is to take the cross product of the vectors formed by the point-segment ([point_x, point_y]-[p_x_k,py_k]) and each segment ([p_x_k+1,py_k+1]-[p_x_k,py_k]).  This will return the direction of the normal vector to the plane formed by the two segments.

As we loop through each segment we check the cross product to see if the angle is within the same range of 0 to pi or pi to 2pi. If one segment is found to be outside the initialialized sign, the point is considered outside the polygon. Let’s look at this loop:

% Check if point lies on same side of each segment.  We will check each
% segment in order.
% Assumption: vertices are given in either clockwise or
% counterclockwise manner.
for k = 1:numNodes
    % determine vectors for each segment around boundary
    % starting between p2-p1, p3-p2, etc.        
    % if at the final node, the vector connects back to p1
    if k == numNodes
        current_poly = [px(1) - px(k), py(1) - py(k), 0];
    else
        current_poly = [px(k+1) - px(k), py(k+1) - py(k), 0];
    end

    %vector from point in space vertex on polygon (point-p(k))
    point_vect = [point(1) - px(k), point(2) - py(k), 0];

    % determine if the point is inside or outside the polygon:
    % if cross product of all polygon_vectors x point_vectors is
    % negative then the point is inside (if polygon is convex) 
    %take cross_product of point vector and polygon vector       
    c = cross(current_poly, point_vect);
    current_sign = sign(c(1,3));
    if k == 1
       check_sign = current_sign;
    elseif check_sign ~= current_sign
        insidePoly = false;
        return
    end   
end  

You’ll notice that for the first segment we need to initialize the sign, with each subsequent segment being compared to this sign. Also, while the vector we are using is 2D, we need to pad the vectors with a 0, in order to perform a cross product, as the cross product provides the perpendicular vector, which requires 3D space.

Here are some test cases which you can try on your own, with a sample graph of the points to be tested given here:

There are four possible conditions: outside, inside, on a polygon segment and on a polygon vertex. Try each for a convex polygon and next time we’ll discuss non-convex polygons.

% Test cases:
% Polygon defined as:
% px = [1,-1,-1,2,3,4];
% py = [-1,-1,1,4,4,1];
% Test point outside polygon. Should return false:
% polygonContainsPoint(px, py, [4,-1])
% Test point within polygon.  Should return true:
% polygonContainsPoint(px, py, [2,2])
% Test point on line segment. Should return false:
% polygonContainsPoint(px, py, [0,-1])
% Test point on vertex.  Should return false:
% polygonContainsPoint(px, py, [1, -1])

If you have any thoughts or questions please don’t hesitate to comment!